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Abstract 

The quality of CMB observations has improved dramatically in the last few years, and will continue to do so in 
the coming decade. Over a wide range of angular scales, the uncertainty due to instrumental noise is now small 
compared to the cosmic variance. One may claim with some justification that we have entered the era of precision 
CMB cosmology. However, some caution is still warranted: The errors due to residual foreground contamination in 
the CMB power spectrum and cosmological parameters remain largely unquantified, and the effect of these errors on 
important cosmological parameters such as the optical depth r and spectral index Ws is not obvious. A major goal for 
current CMB analysis efforts must therefore be to develop methods that allows us to propagate such uncertainties 
from the raw data through to the final products. Here we review a recently proposed method that may be a first step 
towards that goal. 
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1. Introduction 

The great importance of the foreground prob- 
lem for CMB studies has long been recognized 
in the cosmological community, and as a result 
a large number of algorithms have been pro- 
posed, implemented and applied to both simulated 
and real data. Examples are the Maximum En- 
tropy Method llBarreiro et al.ll2004l: iBennett et al. . 
2003bl: iHobson et all llQQSt IStolvarov et all l20o" . 



Internal Linear Combination methods 



llBennett et al.[ bonShl: i Tcgmark et al., 200^^ 
lEriksen et a,1.[l2nn4a»Saha et a,L..2nn,^.Ha,nsen et a,l 



kOOd. iHinshaw et all kmd). Wiener filtering 
llBotichet fc GispertL Il999t IXeemark fc EfstathiouL 
'l996 ) , Independent Component A nalysis meth- 
ods (,Maino et al.. .2002, ■2003: ■DonzcUretall. Em 



2^ 
et all 



and s pectral fitting approaches IjDelabrouille et al.l 
'200^. 

The main emphasis of most of these methods lies 
on establishing an optimal estimate of the CMB sky 
signal, and not so much on quantifying the uncer- 
tainty o n that estimat e . A d ifferent approach was 
taken by iBrandt et al.l l|l994l) who simply adopted 
well-established likelihood parameter estimation 
techniques to solve the problem. Ten years later 
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this work was continued bv lEriksen et alJ l|2006() , by 
taking advantage of recent statistical (in particular 
Markov Chain Monte Carlo) and computational 
developments. 

In this pa per we briefly review the algorithm as 
presented bv Eriksen et al.ll|200(i() . We then describe 
three applications, namely 1) an analysis of simu- 
lated Planck + six-year WMAP data, 2) a prelimi- 
nary analysis of the first-year WMAP data, and 3) 
a preliminary study of experimental design. 



2. Methods 

The most straightforward method for producing 
reliable error bars is provided by standard pa- 
rameter estimation techniques. In particular, the 
Bayesian framework is particularly suited for this 
type of problem. 

Given a set of observed data d and an assumed 
model s{6), 9 being a general set of parameters, 
we simply ask, what is the posterior distribution 
P{9\d)? To answer this, we first recall Bayes' for- 
mula. 



A^c, 



Sii^;e)^ Y, S,{v-6i) 



p{e\d) ex p{d\e)p{e) = c{0)p{0), 



(1) 



where C{9) = P{d\9) is the likelihood and P{0) is a 
prior summarizing our previous knowledge about 9. 
For Gaussian data, the likelihood is given by 

-2 ln£ = x'{0) = [d- s(0)]* C'^ [d - s(0)] , (2) 

where C is the covariance matrix, up to an irrele- 
vant constant. (In the current implementation, we 
assume no pixel-pixel correlations, and C = N is de- 
fined to be the noise covariance matrix.) The prob- 
lem is thus reduced to mapping out the posterior as 
a function of the free parameters 9 using some com- 
putational tool, such as MCMC or grid evaluations. 

This machinery may be applied to microwave 
component reconstruction by noting that different 
signal components have different spectral and spa- 
tial behaviour. By observing the sky in different 
frequencies and directions one may disentangle the 
various contributions and isolate the CMB signal. 

The first step is therefore to choose a suitable 
parametric model. For CMB analysis, a particularly 
convenient choice is that of a sum of independent 
modulated power-law components. 



JVc, 



= J2 A,h{^) — 



(3) 

(4) 



Here one typically would include a CMB term 
(^cmb = Tcnib; fcmh{v) IS the thcrmodynamic-to- 
antenna temperature conversion factor; /3cmb = 0), 
a synchrotron and free-free term [A and/or (3 free 
parameters; /s(j^) = 1), and a thermal dust term. 

These few simple definitions summarize the 
method quite succinctly. No a priori assumptions 
about the sky emission are required, except that 
it can be represented by a parametric model. The 
remaining discussion is concerned mostly about 
how to deal with real-world complications such as 
limited spectral information, low signal-to-noise 
ratio, and computational constraints. It is worth 
emphasizing that even an unrealistically simple sky 
model containing nothing but CMB, free- free, syn- 
chrotron, and dust emission requires at least six 
parameters, and that no CMB experiment to date 
has had even the eight frequencies required to fit 
such a minimal model. Fortunately, the method can 
be extended easily to include other information, as 
will be seen in the next section. 



2.1. Calibration and template fitting 

Two major complications arise when trying to 
apply the above machinery to real- world data. First, 
spectral fits generally require that all sky maps 
are properly calibrated with respect to a common 
zero-point. This is not trivial for any CMB experi- 
ment, and certainly not for differential observatories 
such as the WMAP satellite. Second, as previously 
mentioned, the number of observed frequencies is 
often (i.e., in every experiment performed to date!) 
smaller than the number of parameters one might 
wish to include. For instance, WMAP observes the 
sky at five frequencies, while, ideally, we would like 
to include at least six parameters in a reasonable 
model (CMB, synchrotron, free- free and thermal 
dust amplitudes, and synchrotron and dust spectral 
indicies), still neglecting a possible anomalous dust 
contribution. 

A straightforward way to address both problems 
is to introduce template terms in the signal model. 
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Fig. 1. Marginalized posterior distributions for six parameters fitted to a single pixel by Markov Chain Monte Carlo. This 
example is taken from a simulation corresponding to the combination of Planck and six-year WMAP, for a pixel located exactly 
on the Galactic plane. The vertical lines indicate the true parameter values. (Note that this is not a well-defined quantity for 
the dust spectral index, since we fit for a one-component dust model whereas the simulation is based on a two-component 
dust model.) 



S{i^;9)~>S{n,iy;e) = 

Wpix Wto 



Y^S{i^;e.)+ Y. Mv)g,{i 



(5) 



That is, in addition to the previously defined (and 
usuaUy rather compUcated) spectral models, we also 
fit for a small set of global (and simple) template 
amplitudes. Common templates to include would 
be a monopole term, three dipole terms, and ei- 
ther reasonably well-behaved (e.g,. free-free emis- 
sion through an H^ template) or sub-dominant (e.g., 
thermal dust emission for WMAP) foreground com- 
ponents. For each signal component one is able to 
handle this way, one saves one degree of freedom in 
all subsequent non-linear fits, thus obtaining very 
useful additional stability. 

The cost of this approach is that the problem be- 
comes global, and the computational load is vastly 
increased. At present it is therefore only feasible to 
perform such an analysis at low resolution and with 
a non-linear search. Nevertheless, since the main tar- 
get at this stage is a relatively small number of global 
template amplitudes, little is lost by this approach. 

Note that this extension of the method was dis- 
cussed only briefly by Eriksen et al. ( 2006) in the 
context of monopole and dipole calibration, not gen- 



eral template fits, and it was also not implemented 
at that time. The analysis of the first- year WMAP 
data presented in this paper is therefore the first ap- 
plication of the method. 



2.2. Estimating spectral indices from low-resolution 
maps 

After calibrating and removing well-behaved com- 
ponents from the maps using the above description, 
the next step is to map out the full posterior for 
each pixel with low-resolution maps. The reasons 
for using low-resolution maps at this stage are two- 
fold: On the one hand, estimation of non-linear pa- 
rameters such as spectral indices is quite sensitive 
to noise, and it is therefore highly desirable to have 
high signal-to-noise data. On the other hand, we use 
full-blown MCMC to map out the posteriors for each 
pixel, and with a computational cost of about 100 
CPU-seconds per pixel it is currently not feasible to 
do this for mega-million pixel maps. 

For these two reasons, we smooth each map with a 
wide beam (typically with a Gaussian beam of, say, 
7° FWHM) and downgrade the pixel resolution (typ- 
ically to, say, 2° pixel size; Nsidc = 32 in HEALPix 
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Fig. 2. Simulation constructed to test the algoritlim based 
on a semi-realistic model of the foregrounds: The syn- 
chrotr on emission includes s patial variations in the spectral 
index JGiardino et all 1200211 . and the ther mal dust emission 
follow s the two-component "model 8" of iFinkbeiner et aU 
(19991) . For further information on this simulation, see 
lEriksen etlD i2006D . 

language^ ) prior to the MCMC analysis. The com- 
putational expense is then quite manageable with a 
modern-type cluster, with a total cost of about 500 
CPU hours. 

The final product from this stage is a single joint 
probability distribution for each pixel, and from 
these we find the univariate distributions for each 
parameter by marginalizing over all others. Exam- 
ples of such univariate distributions are shown in 



Figure 1. Finally, we store the posterior mean and 
variance for each parameter and each pixel in the 
form of two sky maps as our final data products. 

2.3. Estimating amplitudes from high-resolution 
maps 

While the procedure described in the previous sec- 
tion gives a full representation of the posteriors for 
the low-resolution sky maps, we also want a good 
representation of the full-resolution posteriors. To 
obtain these, we make the assumption that the spec- 
tral indices vary more slowly than the amplitudes ^ , 
and fix their distributions at the corresponding low- 
resolution values. Given the values of all non-linear 
parameters, the likelihood becomes Gaussian, and 
the analysis may be performed analytically with 
modest compu tational expenses. Fo r details on this 
procedure, see lEriksen et al.l l)2006(l . 

3. Example applications 

In the following sections, we show a few demon- 
strations of how this method performs in prac- 
tice. First, we revie w the simulation described by 
lEriksen et al.l()20n(Tl) . We then show the results from 
a preliminary analysis of the first-year WMAP data. 
We emphasize that these results are preliminary, 
and will be revised with the currently available 
three-year data. Finally, we show some early results 
from a currently on-going study of optimization of 
future experiments. 

3.1. Planck + six-year WMAP simulation 



In order to test the algorithm, lEriksen et alJ 
l|20 06:) constructed a simulation corresponding to a 
combination of six Planck channels (30 to 217 GHz) 
and five WMAP channels (23 to 94 GHz). The 
higher-frequency Planck channels were not included 
because of modelling error confusion associated 
with including multiple dust components. The sim- 
ulation took into account the beam and white noise 
characteristics of each channel separately. Four sig- 
nal components were included: a Gaussian CMB 
realization drawn from the best-fit WMAP ACDM 
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^ Whether this assumption is well justified will only be clear 
after the first high-quality observations (from Planck) are 
analyzed, but based on the limited data available to date, 
this looks to be a reasonable approximation. 
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Fig. 3. Full-sky results obtained from the simulation described in Section 3.1. Shown are reconstructed parameter values (left 
column), actual residual errors (middle column), and estimated uncertainties (right column). The sharp boundaries are caused 
by different assumed models (e.g., free- free is omitted from the fitting model at high latitudes but included at low), which 
leads to order-of-magnitude different uncertainties because of degeneracies. 



power spectrum; a semi-realistic synchrotron com- 
ponent, featuring a sp atially varying spectral index 
l|Giardino et alll2002() : a free-free component with 
fixed spectral index of /3ff = —2.15 I Dickinson et all 



|200S^ : and a two-component thermal dust model 
("model 8" of Finkbeiner et al. 1999). Thus, the 
simulation includes two major complications that 
we must expect to find in real data, namely both 
spatially varying spectral indices and departures 
from simple power law spectra. Three selected 
channel maps are shown in Figure 2. 

These simulations were then analyzed using the 
machinery described above. The results are dis- 
cussed in detail in the original paper; two sample 
results are reprinted here. In Figure 1 we show the 
marginal densities for one single pixel located ex- 
actly on the Galactic plane, and the true values for 



each parameter is marked by a solid line. Clearly, 
the posterior distributions agree very well with the 
input values, both in mean and variance. 

In Figure 3 we show sky maps of the four com- 
ponent amplitudes that were estimated: Posterior 
means are shown in the left column, actual residu- 
als are shown in the middle column, and estimated 
errors are shown in the right column. Again, the re- 
sults agree very well with expectations. 

3.2. The first-year WMAP data 

In this section we present for the first time the 
results from a parametric foreground analysis o f 
the first-year WMAP data l)Bennett et alll2003a|) . 
These data comprise a set of full-sky temperature 
maps at five frequencies (23, 33, 41, 61 and 94 GHz) 



at angular resolutions between 13' and 53' FWHM. 
The noise is assumed to be white, and all beams are 
assumed to be circularly symmetric. 

As described above, the analysis was performed 
in three stages, two at low resolution (7° FWHM 
for WMAP), and the third at high resolution (1°). 
In the first or calibration stage, the free parameters 
were a monopole an d dipole arnplitud e, a free-free 
template amplitude IIFinkbeineil 120031). and a dust 



template amplitude IjFinkbeiner et all 119991) for 
each band, and a CMB temperature and a syn- 
chrotron amplitude and spectral index for each 
pixel. The fit was performed using a non-linear Se- 
quential Quadratic Programming (SQP) algorithm. 

In the second stage, the low-resolution posterior 
distributions are mapped out using MCMC. Only 
CMB and synchrotron (both amplitude and spec- 
tral index) were included in this case. The main re- 
sult from these computations is a synchrotron spec- 
tral index map, shown in the bottom panel of Fig- 
ure 4. Note that the most values lie between —2.0 
and —3.4, which is quite acceptable. (Values larger 
than —2.3 probably indicate residual free-free emis- 
sion, rather than break-down in the synchrotron es- 
timation.) 

In the third stage, we estimate the CMB and syn- 
chrotron amplitudes analytically. The results from 
these calculations are shown in the two panels of 
Figure 4. (The synchrotron amplitude is normalized 
to K band.) Again, the results appear visually quite 
compelling, although a more quantitative analysis 
is warranted. However, we re-emphasize that these 
results are presented only to demonstrate the capa- 
bilities of the method, and not as definitive measure- 
ments. The analysis will soon be revisited based on 
the recently available three-year WMAP data. 

3.3. Experiment design 

The machinery described in Section 2 is also very 
well suited to study optimization of future experi- 
ments. The question we want to answer is, given a 
set of instrumental constraints (e.g., focal plane area 
and power dissipation) how should we distribute the 
number of detectors as a function of frequency in 
order to optimize our ability to extract the CMB 
signal? 

One possible method for answering this question 
would be to apply the above machinery to a wide 
range of allowed instrument configurations, and 
study how the uncertainties (or residuals) change 
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Fig. 4. Results from first-year WMAP analysis. Shown are 
reconstructed CMB signal (top panel), synchrotron ampli- 
tude (middle panel), and synchrotron spectral index (bottom 
panel) . 

with configuration. In a currently on-going project, 
we try to do this in a systematic fashion, and in this 
section we present some very preliminary results 
from this study. 

We consider a generalized version of the simu- 
lations described in Section 3.1, for which the fre- 
quency bands and the noise RMS per frequency may 
be chosen freely. We then impose a set of constraints 
corresponding to the focal plane area of an optical 
system that could fly in space, and feeds at differ- 
ent frequencies of realistic size. Based on these con- 
straints, we generate a grid of possible instrument 
configurations, and run the analysis for a reasonable 
set of these. (Most may be discarded by common 
sense.) 

A few example results are shown in Figure 5: 
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Fig. 5. Example results from the experimental design study. 
The curves represent different configurations allowed by the 
experimental constraints, and show the effective noise per 
frequency channel, which is directly related to the number 
of detectors per frequency. The thickness (and color) indi- 
cates the reconstruction quality: A thicker (and redder) line 
indicates a better reconstruction. The optimal solution here 
is the thick red line, which corresponds to a constant sig- 
nal-to-noise ratio. 



The lines indicate the noise RMS distribution as 
a function of frequency (o'(i^) — crQ{v) / ^/N\iet(y) , 
where iVdct is the number of detectors at frequency 
v) for seven different configurations. The thickness 
(and color) of the lines indicates the resulting un- 
certainty/residual from the component separation: 
A thicker (redder) line means a better reconstruc- 
tion. The optimal distribution of detectors among 
this set is therefore the solid red line. 

From this simple exercise, we may formulate a 
first design principle: The optimal noise distribution 
is the one that corresponds to a constant signal-to- 
noise ratio, where the signal is the sum of CMB and 
foregrounds. Second, the wider the frequency cover- 
age, the better the reconstruction. 

However, there is a major caveat here: This princi- 
ple only holds when any modelling errors (i.e., errors 
due to the fact that the assumed parametric model 
is different from the true signal) are small. Perform- 
ing a similar exercise on simulations with significant 
modelling errors results in a competing principle: 
When the parametric shape of the foreground compo- 
nents are poorly known, the best reconstructions are 
obtained with configurations surrounding the fore- 
ground minimum. The optimal total frequency span 
depends on the magnitude of the modelling errors. 
Both of these principles sound rather obvious, but 
they arc nevertheless worth making explicit. 

In conclusion, two competing effects are at work: 



Reconstruction power is maximized by a wide fre- 
quency coverage, but limited by modelling errors. 
In a future publication, we will quantify these con- 
siderations in much greater detail, and attempt to 
provide some guidelines on the preferred frequency 
ranges for future CMB experiments. Firm conclu- 
sions, however, will require an understanding of 
modelling errors based on better knowledge of the 
foregrounds themselves, knowledge that can come 
only from measurements of the polarized sky over a 
broad range of frequencies with great sensitivity. 



4. Discussion 

As the sensitivity of CMB observations continue 
to improve, the importance of properly character- 
izing the foreground contributions also increases. 
Even with WMAP we are already at a level where in- 
strumental noise is negligible compared to the fore- 
ground ui]£ertaiiittes_^Dvei_a wide range of angular 
scales i Hinshaw et al.L I200(tI) . and this will be even 
more true for Planck and up-coming polarization ex- 
periments. Simple template corrections will be hope- 
lessly inadequate to reach the required sensitivity 
levels. Further, it will be essential to propagate the 
foreground uncertainties through to the final prod- 
ucts, namely the CMB power spectrum and cosmo- 
logical parameters. 

On this background, we argue that the most 
appropriate solution will be based on traditional 
parameter estimation techniques, rather than im- 
age processing techniques. The reason is simply 
that such methods naturally provide the full pos- 
terior distributions, and are integrated easily with 
existing power spectrum and parameter estimation 
methods based on Bayesian parameter estimation. 
An approach that appears particularly promising 
in this respect is that of Gibbs sampling, which 
allows for joint global analysis of foregrounds, 
CMB power spectrum and even cosmological pa- 
rameters (Jewell et al 2004: Wandelt et al. , 200 J: 
iRriksen et a].Ll2n04b[l2noir 



In fact, we end by emphasizing that the particular 
implementation we describe in this paper is of less 
importance than its underlying idea. Many improve- 
ments can be made to the algorithm as such (e.g., 
introducing support for spatial correlation informa- 
tion would be of great value) , but our main conclu- 
sion is independent of such details: CMB component 
separation is a probabilistic problem, and obtaining 
accurate uncertainties is a crucial part of the prob- 



lem. Parameter estimation techniques provide the 
most direct route for doing so. 



References 
Barreiro, R. B., Hobson, M. P., Banday, A. J., 

Lasenby, A. N., Stolyarov, V., Vielva, P., & 

Gorski, K. M. 2004, MNRAS, 351, 515 
Bennett, C. L. et al. 2003a, ApJS, 148, 1 
Bennett, C. L. et al. 2003b, ApJS, 148, 97 
Bouchet, F. R., & Gispert, R. 1999, New Astronomy, 

4, 443 
Brandt, W. N., Lawrence, C. R., Readhead, A. C. S., 

Pakianathan, J. N., & Fiola, T. M. 1994, ApJ, 

424, 1 
Delabrouille, J., Cardoso, J. -P., & Patanchon, G. 

2003, MNRAS, 346, 1089 
Dickinson, C., Davies, R. D., & Davis, R. J. 2003, 

MNRAS, 341 369 
D onzelU, S., et al., 2005, MNRAS, submitted, 

■astro-ph/0507267" 
Eriksen, H. K., Banday, A. J., Gorski, K. M., Lilje, 

P. B. 2004a, ApJ, 612, 633 
Eriksen, H. K., et al. 2004b, ApJS, 155, 227 
Er iksen, H. K., et al. 2006, ApJ, 641, 665, 

|astro- ph/0508268 
Finkbeiner, D. P., Davis, M., & Schlegel, D. J. 1999, 

ApJ, 524, 867 
Finkbeiner, D. P. 2003, ApJS, 146, 407 
Giardino, G., Banday, A. J., Gorski, K. M., Bennett, 

K., Jonas, J. L., & Tauber, J. 2002, A&A, 387, 82 
Gorski, K. M., Hivon, E., Banday, A. J., Wandelt, 

B. D., Hansen, F. K., Reinecke, M., &Bartelmann, 

M. 2005, ApJ, 622, 759 
Hansen, F. K., Banday, A. J., Eriksen, H. K., Gorski, 

K. M., & Lilje, P. B. 2006, ApJ, submitted, 

astro-ph/0603308 
Hinshaw, G., et al. 2006, ApJ, submitted, 

■astro-ph/060345r 
Hobson, M. P., Jones, A. W., Lasenby, A. N., & 

Bouchet, F. R. 1998,MNRAS, 300, 1 
Jewell, J., Levin, S., & Anderson, C. H. 2004, ApJ, 

609, 1 
Maino, D., et al. 2002, MNRAS, 334, 53 
Maino, D., Banday, A. J., Baccigalupi, C., Perrotta, 

F., & Gorski, K. M. 2003, MNRAS, 344, 544 
Sa ha, R., Jain, P ., & Souradeep, T. 2005, 

■astro-ph/0508383" 
Stolyarov, V., Hobson, M. P., Ashdown, M. A. J., & 

Lasenby, A. N. 2002, MNRAS, 336, 97 



Stolyarov, V., Hobson, M. P., Lasenby, A. N., & Bar- 
reiro, R. B. 2005, MNRAS, 357, 145 

Tegmark, M., & Efstathiou, G. 1996, MNRAS, 281, 
1297 

Tegmark, M., de Oliveira-Costa, A., & Hamilton, 
A. J. 2003, Phys. Rev. D, 68, 123523 

Wandelt, B. D., Larson, D. L., & Lakshmi- 
narayanan, A. 2004, Phys. Rev. D, 70, 083511 



